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1 Introduction 


This report describes the mathematical model and associated code used to simulate a high speed 
civil transport aircraft - the Boeing Reference H configuration. The simulation was constructed to 
support advanced control law research by providing an appropriate tool for control system 
engineering. For this reason, the simulation was implemented in a commercially available and 

widely used controls engineering software package called M ATL AB/SIMULINK® . 

To date there have been 3 formal releases of the Reference H simulation math model as shown in 
the table below: 

Cycle 1 September, 1994 (Dawdy et al- June, 1994) 

Cycle 2 A March, 1995 (Stephens et al - March, 1995) 

Cycle 2B July, 1995 (Stephens et al - July, 1995) 

This report describes SIMULINK implementation of the Cycle 1 Reference H simulation model. 
Future versions of this report will be issued on a regular basis to document subsequent releases of 
the simulation math model. 

MATLAB and its toolboxes provide the linear algebra, matrix computation, numerical analysis 
capabilities, and advanced control law analysis and synthesis methods that are suitable for control 
system engineering. SIMULINK provides a graphical, block diagram format for constructing and 
simulating dynamic systems. The simulation software developed on the MATLAB/SIMULINK 
platform can be used for obtaining nonlinear dynamic responses of the aircraft and also provides 
the capability to find trim solutions and linear models at specified flight conditions. Thus, the 
developed simulation provides the necessary basis for developing and testing candidate flight 
control systems for a representative high speed civil transport aircraft. 

The simulation math model includes the nonlinear, six degree-of-freedom equations that govern the 
motion of a rigid aircraft in atmospheric flight over a rotating spherical Earth. A Dryden turbulence 
model and the 1962 Standard Atmosphere model are provided to simulate the Earth atmosphere. 
The Reference H configuration is modeled using the Boeing cycle 1 simulation data base and 
consists of a tabular aerodynamic model, an engine model, and a mass model. The aerodynamic 
model has the option to include the effects of quasi-static elastic deformation on the six rigid-body 
aerodynamic forces and moments. Landing gear, flight control and the actuator models are not 
documented in this report will be covered in future releases of the simulation and itsdocumentation. 

The equations of motion used in the simulation are presented in Chapter 2. Chapters 3, 4, and 5 
present brief descriptions of the Reference H configuration aerodynamic, engine, and mass models 
used in the simulation. The equations used to compute the accelerations at a given sensor location 
are presented in Chapter 6. The atmosphere model and the turbulence model implemented in the 
simulation are described in Chapter 7. Chapter 8 describes the overall simulation layout and the 
procedure to run the nonlinear simulation with an example case. The trimming procedure and linear 
model generation procedures with examples are presented in Chapters 9 and 10, respectively. 

To validate the developed software, a dynamic check case was generated using the MATLAB/ 
SIMULINK based simulation and compared to dynamic responses with those generated from 
independent simulations implemented by Boeing Commercial Aircraft Company and by the real- 
time simulation facility support group at NASA Langley Research Center. The results of the 
comparison are also included. 


® MATLAB and SIMULINK are registered trademarks of The MathWorks, Inc. 
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2 Aircraft Equations of Motion 


2.1 Underlying Assumptions 

The equations governing the motion of a rigid aircraft in atmospheric flight are available from 
several sources (Etkin - 1972 and McFarland - 1975), and are not derived in this paper. The 
equations discussed in this paper are taken from Etkin - 1972. The equations utilize a body 
reference frame with the origin at the aircraft center of gravity, schematically shown in Figure 2. 1 . 
The positive x axis passes through the aircraft nose. The positive y axis passes through the right 
wing and the positive z axis passes through the bottom of the aircraft fuselage. The equations of 
motion were developed based on the following assumptions: 

1 . The Earth is a sphere rotating with constant rotation (0 E on an axis fixed in inertial 
space. 

2 . Acceleration due to gravity is a radial vector which acts at the aircraft mass center. 

3. The aircraft is a rigid body which is symmetrical about the plane formed by the body 
frame x axis and z axis. 


2.2 Translational Equations of Motion 

Let the translational inertial velocity of an aircraft, expressed in body frame components, be given 
by [ u v w]' and let the inertial rotational velocity of an aircraft in body frame components be given 
by [p q r]'. The rotational velocity of the Earth in body frame components is given by 


Pb 


cosA 

<lb 

” - [Tflv} 

0 

kJ 


—sin A 


( 2 . 1 ) 


where A is aircraft latitude and [r sv -] is the transformation matrix from a vehicle-carried local 
horizontal reference frame to a body reference frame, which is given by 


cos# cos yr 

cos#sinyr 

-sin# 

sin 0 sin# cos yr 

sin <p sin# sin yr 

sin <j> cos# 

— cos^sinyr 

+ cosdcosyr 


cos 0 sin# cos yr 

cos <p sin# sin yr 

cos 0 cos# 

+sin0sinyr 

-sin^cosyr 




where (p , #, and yr are the Euler roll angle, pitch angle, and yaw angle respectively. The 
vehicle-carried local horizontal reference frame is attached to the aircraft with the origin at the center 
of gravity. The positive x axis points to the north, the positive y axis points to the east, and the 
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positive z axis points vertically downward toward the Earth’s center. Let m be the mass of the 
aircraft and let [ F x F y F z }' represent the total external force components acting on the aircraft. The 
equations describing translational motion are given by 



u + (q + <lb)w-(r + r£)v 


F x 


v + (r + rb)u-(p + Pt)w 

> = < 

Fy • 


w + (p + Pb)v-(q + qb)u 


Fz. 


where the total force vector is given by 


F x 

Fy 

► — < 

F x 

Fy 

> + « 


> + ^ 

F x 

Fy 

Fz. 


Fz. 

Aero 

F z 

Prop 

fz. 


(2.3) 


(2.4) 


The subscripts ‘Aero’, ‘Prop', and ‘ Grav ’ denote aerodynamic, propulsive, and gravitational 
forces respectively. The components of the gravitational force vector are 


(2.5) 


[F 1 

r x 


-sinO 

Fy 

> = mg- 

cos^sin^ ► 

Fz. 

Grav 

cos 8 COS 0 


Gravitational acceleration is denoted by g and is given by 

s = *„(« £ /R) 2 . 


( 2 . 6 ) 


where g Q is the gravitational acceleration at sea level, R E is the radius of the Earth, and R is the 
distance from the center of the Earth to the aircraft mass center. Equation (2.3) can be rewritten as 
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0 
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F x 
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(2.7) 


2.3 Rotational Equations of Motion 

Let the total external moments, in body frame components, acting on the aircraft be given by 
[ M x M y M z ]'. I x is the roll moment of inertia, I y is the pitch moment of inertia, I z is the yaw 
moment of inertia and I K is the product of inertia. The equations describing the rotational motion 
of an aircraft are given by 
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( 2 . 8 ) 


Ph - (<1P + r)lxz + qr[l z - I y ) 


M x 

ql y +(p 2 -r 2 )l xz + pr(l x -I z ) 
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My ► 

rI z +{qr-p)I xz + pq(l y -I x ) 


M z. 


where the total moment components are given by 
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M x 

My 

► = < 

My 

> + 4 

My > 

M z 


Mz. 

Aero 

Mz. 


(2.9) 


Equation (2.8) neglects the rotor terms and the centripetal acceleration associated with the Earth’s 
rotation. The time rate of change of the moments of inertia are assumed to be negligible. 
However, the moments of inertia can be varied parametrically slowly over time to reflect fuel 
bum. Equation (2.8) can be rewritten as 
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2.4 Combined Accelerations 

When the body frame inertial velocity components are combined to form the generalized velocity 
vector given by 


u 

V 



' 7 


IM 

equation (2.7) and equation (2.10) can be combined to form 

[M],{V}+[QIM],{V} = {F}. 
The inertial mass matrix [M] ; is given by 


( 2 . 11 ) 


( 2 . 12 ) 
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m 0 0 0 0 0 

0 77i 0 0 0 0 

0 0 772 0 0 0 

[M]l= 0 0 0 I x 0 -/„ 

0 0 0 0 I y 0 

0 0 0 -/« 0 I z 

and [Q] represents the body frame inertial rotation, which is given by 

0 -(r+r/) (q + qf) 0 0 0 

(r+r/) 0 ~{ p + p !) 0 0 0 

[£il = + (P + rf) 0 0 0 0 

1 0 0 0 0 -r q 

0 0 0 r 0 -p 

0 0 0 — q p 0 

and {i 7 } is the generalized force vector, given by 

V 

Fy 

F Z 

m/' 

My 

M Zj 



(2.13) 


(2.14) 


(2.15) 


2.5 Explicit Acceleration Solution 

Like most aerodynamic models, the Reference H aerodynamic model generates forces and 
moments which contain unsteady flow effects. The unsteady flow contributions to the 
aerodynamic forces and moments are functions of acceleration. Equation (2.12) must be modified 
in order to obtain an explicit definition for the aircraft acceleration. Therefore, the force and 
moment contributions from unsteady flow must be removed from the generalized force vector 
defined in equation (2.15). This removal is accomplished by separating the aerodynamic forces 
and moments according to whether the flow is steady or unsteady. To see how this is done, first 
consider the longitudinal force and moment coefficients generated by the Reference H aerodynamic 
model. These force and moment coefficients contain linear contributions from non- 
dimensionalized forms of the time rate of change of angle of attack and the time rate of change of 
pitch rate. Defining ri as the time rate of change of angle of attack and q as the time rate of change 
of pitch rate, the longitudinal force and moment coefficients, in stability frame components, are 
expressed as 
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(2.16) 


Cy — ... + Cy . 
X Xa 


f ac ^ 


K^TJ 


+c xA 


f ~ \ 2 
C 

<^t) ’ 


where ‘ %' can be replaced by ‘ £>’ for drag, * U for lift, or ‘ m ’ for pitching moment. The mean 
aerodynamic cord is given by c and the aircraft speed relative to the air mass is given by V T . 

Similar to the longitudinal case, the lateral force and moment coefficients depend linearly on non- 
dimensionalized forms of the time rate of change of sideslip angle, the time rate of change of roll 

rate and the time rate of change of yaw rate. Let )3 denote the time rate of change of sideslip angle, 
p denote the time rate of change of roll rate, and r denote the time rate of change of yaw rate. 

The lateral force and moment coefficients, in stability frame, are given by 




+ 


(jb} 

\ 2V TJ 


\ + C. 


x-P I 


\ 2 


y2V T j 


+ C Xf n 




2V, 




(2.17) 


where, for the lateral case, ‘ %' is replaced by ‘ Y' for sideforce, ‘ /’ for rolling moment, or ‘ n’ 
for yawing moment. The aircraft wing span is denoted by b . The steady and unsteady flow 
aerodynamic forces and moments are transformed from the stability frame to the body frame before 
being added to the generalized force vector. This transformation process is described in Chapter 3. 

To begin separating the aerodynamic forces and moments, first consider the angle of attack in still 
air, which is defined as 


a - tan 



(2.18) 


The time derivative of equation (2.18) is 


a = 


uw-wu 


The aircraft sideslip angle in still air is given by 


P 


- sin 


-l 




\ V TJ 


(2.19) 


(2.20) 


The time derivative of equation (2.20) is 


p_ ^v-vVr 
Vj-Vm 2 + w 2 


( 2 . 21 ) 


Substituting equations (2.19) and (2.21) into equations (2.16) and (2.17) respectively, the 
aerodynamic forces and moments are clearly seen to be linear functions of body frame 
accelerations. Now, the generalized aerodynamic force vector {F} Aero can be expanded using 
Taylor Series as follows: 
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{f h«„ = w 


<?{F} 


Aero 


Aero 


{*}-» 3{v\ 


{ l >]-0 


w- 


1 <? 2 {f) 


Aero 


21 


m 


M 2 


(2.22) 


{v}-0 


The linear dependence of the aerodynamic forces and moments on body frame accelerations 
simplifies equation (2.22) to 


{n^r, = {n 


3{F} 


Aero 


Aero 


W=° d{v} 


M=° 


M 


(2.23) 


The generalized aerodynamic force vector is now separated into steady and unsteady components. 
Equation (2. 12) can now be written as 

W],{v} = +{F} 

Grav jyJ = Q -[M] A {y}. (2.24) 


where 







(2.25) 


[M] a represents an increment to the inertial mass matrix [M] / , and {F} Crav consists of the gravity 

force vector defined in equation (2.5) and zero moment contributions to make the generalized 
gravity vector. Defining the generalized inertial force vector as 

{F), = (2-26) 

equation (2.24) can be written as 

(M* +M,)M = + {F} Prop + {F} Crm + {F} „ J{v}, 0 - (2.27) 

Finally, let the total mass matrix [M] be given by 


[M] = [M] a +[M],. (2.28) 

The equations of motion as implemented in the simulation are given by 

[M]{V] = {F},+ {F} en>p + {F} Crav + {F} A (2.29) 

At each integration time step equation (2.29) is solved for the unknown body frame accelerations 
by left-multiplying each side of the equation by the inverse of the total mass matrix. 
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2.6 Aircraft Mass Center Position 

Let the aircraft mass center position relative to a round and rotating Earth be given by [ T A h ] ' 
where T defines longitude, A defines latitude, and h defines altitude above sea level. Also, let 
the aircraft inertial translational velocity in local frame components be given by [ V N V E V D ]'. The 
local horizontal reference frame velocity components are given by 


V N 

k - \Tbv\ ' 

r ^ 

U 

<V E 

V 

V D. 


w 


The time rate of change of the aircraft position vector is given by 


(2.30) 


T 


VeJRc osA' 

< A 


v n /r > 

h 

' 

. -V D , 


(2.31) 


where R is the distance from the center of the earth to the aircraft mass center and V Ec is the 
eastward velocity component corrected for Earth rotation, which is given by 


V E c = V E - (O e Rc os A. 


(2.32) 


A singularity occurs in equation (2.31) when A is equal to 90°. The calculation of cos A in 
equation (2.31) is therefore implemented in the simulation as 


cosA = 


A 

£ w 

cos A 


cosA| < £ 
cosA| > e 


, where £ = 10 


-10 


(2.33) 


Division by zero in equation (2.3 1) is therefore avoided. Equation (2.31) is integrated at each 
simulation time step to yield the aircraft position vector. 


2.7 Aircraft Orientation Using Euler Angles 


The orientation of the aircraft body reference frame with respect to a vehicle-carried local horizontal 
reference frame is described by the Euler angles <j>, 0, and if/. Let the relative rotational velocity 
between the body reference frame and the local horizontal reference frame be given by [ P Q R ] ' . 
The time derivatives of the Euler angles are then calculated as 



P + 


1 

COS0 


( Q sin 0 sin 0 + R cos 0 sin 0) 
Qcostp- /?sin0 


(<2sin 0 + Rcos<j>) 

COS0 


(2.34) 
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where 


p 


r ^ 

p 


(<o £ + t)cosA 

Q 

* 4 

q 

-[T.A- 

-X 

R 

K J 


r 


-{(0 E + r)sinA 


(2.35) 


A singularity occurs in equation (2.34) when 6 is equal to 90°. The calculation of cos0 in 
equation (2.34) is implemented in the simulation as 

|cos0| < £ ,» 

1 1 , where e = 10 °. (: 

|cos0| > e 

Division by zero in equation (2.34) is therefore avoided. Equation (2.30) is integrated at each 
simulation time step to yield the Euler angles. 

2.8 Flight Path Calculations 


cost? = < 


e 

£ |0| 

cost? 


Horizontal flight path angle, y h , and vertical flight path angle, y v , are given by 


Yh ~ tan 


-if V, 


V N J 


(2.37) 


y v = tan 


-i 


-V„ 


.Jvl+vZ 


(2.38) 


2.9 Summary 

Let the simulation state vector be defined as x = \u, v, w, p, q, r, h, X, T, (p, 6, y/] T . The 

equations required to propagate the state vector from t to t + At are summarized in this section. 
The translational and rotational equations of motion of the aircraft are given by equation (2.29). 

The total mass matrix in equation (2.29), [M], has two components, the inertial mass matrix [M] / 

and the apparent mass matrix [M] a . The inertial mass matrix [M] { is defined by equation (2.13) 
and the computation of the apparent mass matrix [M] A is described in section 3.5. The inertial 
generalized force vector, {F} ; , is defined in equation (2.26) and the gravity force vector by 
equation (2.5). The calculation of the generalized aerodynamic force vector, {F} Ae J{v} =0 * i s 

described in section 3.4 . The computation of generalized propulsion force vector, {F } Prop , is 

described in Chapter 4. The time rate of change of the aircraft’s center of mass position is given by 
equation (2.31). The state vector and the transformations given by equations (2.30), (2.32), and 
(2.2) define the necessary quantities to propagate the aircraft’s position. The time rate of change of 
the aircraft’s orientation with respect to the local horizontal reference frame is given by equation 
(2.34). The state vector, equations (2.35), (2.31), and (2.2) define the necessary quantities to 
propagate the aircraft’s angular position. 
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3 Aerodynamic Model 


3.1 Aerodynamic Model Description 

The Boeing Reference H cycle 1 aerodynamic force and moment coefficient buildup is discussed in 
detail in Dawdy et al - 1994. The rigid airframe aerodynamic force and moment coefficients are 
constructed from a collection of intermediate terms which are calculated from interpolated lookups 
of tabular data. The aerodynamic model includes a quasi-static aeroelastic logical variable, 

RIGID. If RIGID is set to false then quasi-static aeroelastic adjustments for flexibility are included 
in the aerodynamic force and moment coefficient buildup. 

The aerodynamic model implemented in the simulation is divided into two parts. The steady flow 
part performs the aerodynamic force and moment coefficient buildup excluding the effects due to 
the rotational acceleration components, the time rate of change of the angle of attack, and the time 
rate of change of the sideslip angle. The calculation of steady flow aerodynamics is desribed in 
section 3.4. The unsteady flow part of the aerodynamic model accounts for the effects rotational 
acceleration, the time rate of change of the angle of attack, and the time rate of change of the 
sideslip angle with the apparent mass matrix and is desribed in section 3.5. 

All aerodynamic terms are resolved about the stability axis and are referenced to the aerodynamic 
reference center which is defined in Table 3.1. 

Table 3. 1 Aerodynamic Reference Center Coordinates 
Fuselage Station (+ aft) 2106. in 

Wing Buttock Line (+ right) 0. in 

Water Line (+ up) 226. in 


3.2 Control Surfaces 

Figure 3.1 defines the control surfaces for the cycle 1 aerodynamic model. The aerodynamic 
effects of the spoilers; control surfaces DSP 1, DSP2, DSP3, and DSP4 shown in Figure 3.1, 
are not included in the cycle 1 aerodynamic model. The control surface deflection sign convention 
is given in Table 3.2. 


3.3 Aerodynamic Force and Moment Calculations 

The equations of motion described in Chapter 2 are derived for forces and moments acting at the 
aircraft center of gravity. As was mentioned previously, the cycle 1 aerodynamic model calculates 
forces and moments acting at the aerodynamic reference center. Following the work done by 
Buttrill et al - 1992, the inputs to the aerodynamic model are first transferred from the center of 
gravity to the aerodynamic reference center. The aerodynamic model is called and the aerodynamic 
forces and moments are then transferred from the aerodynamic reference center back to the center 
of gravity. The process of transferring the aerodynamic model inputs from the center of gravity to 
the aerodynamic reference center is described below. 
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Figure 3.1 Reference H Configuration 




Table 3.2 Control Surface Sign 

Convention 


Control Surface Description 

Positive Deflection 



left outboard leading edge flap 

leading edge down 


DLE2 

left inboard leading edge flap 

leading edge down 


DLE3 

right inboard leading edge flap 

leading edge down 


DLE4 

right outboard leading edge flap 

leading edge down 


DTE1 

left outboard flaperon 

trailing edge down 


DTE2 

left mid flaperon 

trailing edge down 


DTE3 

left inboard flaperon 

trailing edge down 


DTE4 

left root flap 

trailing edge down 


DTE5 

right root flap 

trailing edge down 

- 

DTE6 

right inboard flaperon 

trailing edge down 


DTE7 

right mid flaperon 

trailing edge down 


DTE8 

right outboard flaperon 

trailing edge down 


DRUD1 

lower rudder segment 

trailing edge left 


DRUD2 

mid rudder segment 

trailing edge left 


DRUD3 

upper rudder segment 

trailing edge left 


DSTAB 

stabilizer 

trailing edge down 


DELEV 

elevator 

trailing edge down 


DVF 

vortex fence 

trailing edge up 



First the air relative translational velocity at the center of gravity is calculated by subtracting the 
velocity of the air mass due to wind, gusts, and turbulence from the inertial translational velocity 
of the aircraft at the center of gravity. The relative velocity components at the center of gravity are 
given by 


V 


r > 

u 


r y 

u airmass 


> = < 

V 

> — i 

v • > 

v airmass 

."a. 


w 

V. J 


y^airmass. 


The relative translational velocity components are transferred to the aerodynamic reference center 

by 


U arf 


u a 


' 0 

— r 

q " 

x arf 

v arf 

> = < 

v a 

> + 

r 

0 

-p 

yarf > 

w arf 





P 

0 

^arf ^ 


where [ x ar f y ar f z ar f]' defines the distances from the aircraft center of gravity to the aerodynamic 
reference center along body frame axes. Positive x ar j locates the aerodynamic reference center 
ahead of the center of gravity, positive y ar j locates the aerodynamic reference center to the right of 
the center of gravity and positive z ar f locates the aerodynamic reference center below the center of 
gravity. The angle of attack at the aerodynamic center is computed as follows: 
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(3.3) 


°v = tan 


-1 


fw a^ 


\ U arf ) 

The sideslip angle at the aerodynamic reference center is given by 

( \ 


Parf= S^ 


-1 


Vf 

Vt 

V T mf J 


where 


(3.4) 



+ v; 


arf 


+ W. 


arf 


The Mach number at the aerodynamic reference center is given by 



(3.5) 


(3.6) 


where a is the speed of sound. Dynamic pressure is left referenced to the center of gravity and is 
calculated as 


i = (3.7) 

where p is the air density in slugs/ft 3 and V Ta is relative airspeed at the center of gravity. The 

body frame rotational velocities at the aerodynamic reference center are assumed to be equal to the 
rotational velocities at the center of gravity. 


3.4 Steady Flow Outputs 


The quantities defined in equations (3.3) through (3.7) are combined with the altitude of the center 
of gravity, aircraft gross weight, body frame rotational velocities at the aerodynamic reference 
center, and control surface deflections to form the inputs to the aerodynamic model. The 
aerodynamic model defines the steady flow aerodynamic force and moment coefficients which are 
then used to calculate the following stability frame force and moment components: 
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(3.8) 
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where S defines the aircraft wing reference area, the symbol * SA ’ denotes stability frame, and the 
symbol ‘ sf ' denotes steady flow. The stability frame components of the forces are then converted 
to body frame components as follows: 


F x 

BA 

F x 

F ', 

ii 

V 

> 

arf 

kJ 


aif 


(3.9) 


where the symbol ‘ BA' denotes body frame and [T gs ] is the transformation matrix from a stability 
reference frame to a body reference frame and is given by 



cos a ar} 
0 


sina^ 


0 -sina^ 

1 0 

0 cos a arf 


(3.10) 


Similarly, the body frame components of the steady flow aerodynamic moments are given by 


M x 

My 

M 7 


BA 


=fr»} 


arf 


M x 

My 

My 


SA 


arf 


(3.11) 


The body frame components of the steady flow aerodynamic forces at the center of gravity are 
equal to those at the aerodynamic reference center. The body frame components of the steady flow 
moments at the aircraft center of gravity are given as 
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(3.12) 


Finally, the generalized steady flow aerodynamic force vector discussed in Section 2.5 is given by 


i F } A ero 




(3.13) 
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3.5 Unsteady Flow Outputs 

The aerodynamic effects of unsteady flow axe included in the apparent mass matrix [M ] A . This 
section describes the calculations necessary to define the elements of the 6x6 apparent mass matrix. 
Taking the derivative with respect to time of equation (3.2) yields 


a s ' 


U a ~ Warf ^ qZarf 


> = - 

Va+rXarf- PZarf > 

r- 

** 


w a -qx arf + py a rf 


(3.14) 


The time rate of change of the angle of attack at the aerodynamic reference center is given by 

u arf™arf ~ w arf“arf 


Uarf = 


u arf W arf 


(3.15) 


and the time rate of change of the sideslip angle at the aerodynamic reference center is given by 

^T^arf ~ v arf^T 


Parf ~ ' 


-a/m 2 + w 2 


(3.16) 


Substituting the expressions for u at f , v aT j, and w ar j given by equation (3.14) into equation 
(3.15) and equation (3.16) yields 


V-arf = 


_ u arf{< - 4*arf + py a rf)- w arf( il a ~ Warf + Warf) 


u arf + w arf 


(3.17) 


and 


Parf ~ 


V a + rx arf - pz arf V ar f(u a rfU arf + VgrfVarf + WgrfWgrf) 

-Ju, 


arf ^arf 


r2 f? 


(3.18) 


T V u arf " r w arf 


Transforming the a and q stability derivatives from stability frame components to body frame 
components yields 


Cx n - $ ma arf - C Dn cosa arf 


Cz n = -C^ COS cc arf - C Dri sin a arf . 


(3.19) 

(3.20) 


where * T] ’ is replaced by 4 a ’ or ‘ q ’ . The body frame components of the unsteady flow 
aerodynamic forces at the aerodynamic reference center are given below: 


Fy j — qS Cy . 

*usf ^ A <*arf 


f ^arf^ 
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+ qSC x .q 
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V 2V t. 


(3.21) 
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where the symbol ‘ usf ’ denotes unsteady flow. The expressions for a ar j and fi ar f given in 
equation (3.17) and equation (3.18), respectively, can be substituted into the equations (3.21) 
through (3.23) and equations (3.26) through (3.28) to yield the unsteady flow force and moment 
components as functions of body frame accelerations. Transferring the unsteady flow forces from 
the aerodynamic reference center to the center of gravity gives 



The unsteady flow aerodynamic moments in body frame components at the center of gravity are 
given by 
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The apparent mass matrix, defined by equation (2.25), is represented as 
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(3.31) 


where the analytical expressions for the elements of the apparent mass matrix are defined in 
Appendix A. 


3.6 Aerodynamic Model Implementation 

The aerodynamic model is a FORTRAN S-fimction which is compiled as a MATLAB FMEX file. 
The aerodynamic model FMEX file is dynamically linked into the simulation during run time. The 
input parameters for the aerodynamic model are given in Table 3.3. INIT_FLAG is an 
initialization variable. If INIT_FLAG is set equal to 1, all the cycle 1 aerodynamic data is read 
into memory. If INTr_FLAG is 0, then the data is not read into memory. The variable 
RIGID_FLAG is used to set the value of the logical variable RIGID. Setting RIGID_FLAG equal 
to 1 assigns RIGID to true. If RIGID_FLAG is set equal to 0, RIGID is false. Table 3.4 lists the 
aerodynamic model outputs. 
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Table 3.3 Aerodynamic Model Input Parameters 


Variable 

Description 

Symbol 

INIT FLAG 

initialization flag 


RIGIDJFLAG 

quasi-static aeroelastic flag 


U_ARF 

x-body component of inertial translational velocity at arc (fit/s) 

u arf 

V_ARF 

y-body component of inertial translational velocity at arc (fit/s) 

V arf 

W_ARF 

z-body component of inertial translational velocity at arc (ft/s) 

w arf 

PBA_AR 

x-body component of rotational velocity at arc (ft/s) 

P 

QBA_AR 

y-body component of rotational velocity at arc (ft/s) 

q 

RBA_AR 

z-body component of rotational velocity at arc (ft/s) 

r 

X_ARF 

distance from eg to arc along x-body axis (ft) 

x arf 

Y_ARF 

distance from eg to arc along y-body axis (ft) 

y<irf 

Z_ARF 

distance from eg to arc along z-body axis (ft) 

Zarf 

ALPDEG 

angle-of-attack at arc (degrees) 

°w 

BETADEG 

sideslip angle at arc (degrees) 

Parf 

VTOTAL 

airspeed at the arc (ft/s) 

V T 

MACH 

Mach number 

Marf 

QBAR 

dynamic pressure at eg (lbf/ftO 

<7 

ALT 

altitude of eg (ft) 

h 

DLE1 

left outboard leading edge flap deflection (degrees) 

DLE1 

DLE2 

left inboard leading edge flap deflection (degrees) 

DLE2 

DLE3 

right inboard leading edge flap deflection (degrees) 

DLE3 

DLE4 

right outboard leading edge flap deflection (degrees) 

DLE4 

DTE1 

left outboard flaperon deflection (degrees) 

DTE1 

DTE2 

left mid flaperon deflection (degrees) 

DTE2 

DTE3 

left inboard flaperon deflection (degrees) 

DTE3 

DTE4 

left root flap deflection (degrees) 

DTE4 

DTE5 

right root flap deflection (degrees) 

DTE5 

DTE6 

right inboard flaperon deflection (degrees) 

DTE6 

DTE7 

right mid flaperon deflection (degrees) 

DTE7 

DTE8 

right outboard flaperon deflection (degrees) 

DTE8 

DRUD1 

lower rudder segment deflection (degrees) 

DRUD1 

DRUD2 

mid rudder segment deflection (degrees) 

DRUD2 

DRUD3 

upper rudder segment deflection (degrees) 

DRUD3 

DSTAB 

stabilizer deflection (degrees) 

DSTAB 

DELEV 

elevator deflection (degrees) 

DELEV 

DVF 

vortex fence deflection (degrees) 

DVF 

DGEAR 

landing gear extension (degrees) 

DGEAR 

GW 

gross weight (lbf) 
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Table 3.4 Aerodynamic Model Output Parameters 


Variable 

Definition 

Symbol 

MA 

apparent mass matrix (see Appendix A) 

Ma 

FXBA_FM 

x-body component of steady flow aerodynamic forces (lbf) 

% 

FYBAJFM 

y-body component of steady flow aerodynamic forces (lbf) 

Fy 

V 

FZBA_FM 

z-body component of steady flow aerodynamic forces (ft-lbf) 

F z 

LMBAJFM 

x-body component of steady flow aerodynamic moments (ft-lbf) 


MMBA_FM 

y-body component of steady flow aerodynamic moments (ft-lbf) 

M Y„ 

NMBA_FM 

PSA_AR 

QSA_AR 

RSA_AR 

z-body component of steady flow aerodynamic moments (ft-lbf) 
x-stability component of rotational velocity at arc (ft/s) 
y-stability component of rotational velocity at arc (ft/s) 
z-stability component of rotational velocity at arc (ft/s) 

M z* 

CD 

steady flow drag coefficient at arc 


CL 

steady flow lift coefficient at arc 

C L 

L sf 

CM 

steady flow pitching moment coefficient at arc 

C -V 

CR 

steady flow rolling moment coefficient at arc 

c v 

CY 

steady flow sideforce coefficient at arc 

Cy f 

CN 

steady flow yawing moment coefficient at arc 
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4 Engine Model 


4.1 Engine Model Description 

The engine model used in the simulation is representative of a supersonic transport engine. The 
engine is a mixed flow turbofan which incorporates a downstream mixer, 50 % aspiration ratio, a 
nozzle, and a mixed-compression translating centerbody inlet. All four engines are shown in 
Figure 3.1. The model data covers altitudes from sea level to 70,000 feet and Mach numbers from 
0 to 2.4. Gross thrust and ram drag are not included in the model; rather, installed net thrust is 
used. Inlet forces are not specifically modeled, therefore this model does not realistically simulate 
engine failures. The engine model is discussed in Dawdy et al - 1994. 


4.2 Engine Model Dynamics 

The engine model simulates second order nonlinear engine dynamics. The engine dynamics are 
shown in Figure 4.1. The cockpit power lever angle range is limited to values between 0° and 

100° . The limited power lever angle is then directed into a lead-lag filter followed by a rate limiter 
and a lag filter. The filtered cockpit power lever angle value from each throttle is then used to 
obtain the net installed thrust for each engine. 


4.3 Engine Model Buildup Equations 

The engine model calculates the steady state or instantaneous installed net thrust produced by each 
engine for a given power setting. Maximum and minimum installed net thrust limits are calculated 
from interpolated tabular data which is a function of altitude and Mach number. The thrust limits 
are then used with filtered cockpit power lever angle to compute installed net thrust per engine. 
The total installed net thrust from all four engines is used to compute the total thrust-induced body 
frame forces and moments about the aerodynamic reference center. The aerodynamic reference 
center location is defined in Table 3.1. 


4.4 Engine Model Implementation 

The engine model is implemented in two parts. The transfer function implementation of the engine 
dynamics is shown in Figure 4.2. Cockpit power lever angle for each throttle is fed into the engine 
dynamics. The power lever angle signals are biased to zero using the initial cockpit power lever 
angle values. The bias is used to initially set the signal to zero since the filter states are initially set 
to zero in SIMULINK. The biased cockpit power lever angle signals are limited accordingly and 
then directed into the lead-lag filter, followed by the rate limiter, and the lag filter. The bias is 
then removed from the signals, yielding the filtered cockpit power lever angles which are then 
input to the second part of the engine model. The second part of the engine model calculates 
thrust-induced body frame forces and moments about the aerodynamic reference center. This part 
of the engine model is a FORTRAN S-function which is compiled as a MATLAB FMEX file. The 
FMEX file is dynamically linked to the simulation during run time. The engine model input 
parameters are given in Table 4. 1 . The variable TLIMIT_FLAG is used to set the logical variable 
TLIMIT. If TLIMrr_FLAG is set equal to 1, TLIMIT is set to true. Setting TLIMTT to true will 
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Saturation Rate Limit 
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Figure 4.2 Reference H Cycle 1 Engine Dynamics SIMULINK Implementation 






cause the engine model to calculate maximum and minimum installed net thrust as a function of 
altitude and Mach number. If TLIMIT_FLAG is set equal to 0, then TLIMIT is set to false and 
maximum installed net thrust is set equal to 100,000 lbf and minimum installed net thrust is set to 0 
lbf. Table 4.2 gives the engine model output parameters. 


Table 4. 1 Engine Model Input Parameters 

Variable 

Description 


Symbol 

Part 1 PLAl.CP 
PLA2.CP 
PLA3.CP 
PLA4_CP 

throttle 1 cockpit power lever angle (deg) 
throttle 2 cockpit power lever angle (deg) 
throttle 3 cockpit power lever angle (deg) 
throttle 4 cockpit power lever angle (deg) 


PLAl.CP 

PLA2_CP 

PLA3_CP 

PLA4_CP 

Part 2 TLIMTTJLAG 
HPRESS_AT 
MACH.AR 
PLA1F_EN 
PLA2F_EN 
PLA3F_EN 
PLA4F_EN 

maximum and minimum installed net thrust flag 
altitude (ft) 

Mach number 

throttle 1 filtered cockpit power lever angle (deg) 
throttle 2 filtered cockpit power lever angle (deg) 
throttle 3 filtered cockpit power lever angle (deg) 
throttle 4 filtered cockpit power lever angle (deg) 

h 

M 

PLA1F_EN 

PLA2F_EN 

PLA3F_EN 

PLA4F_EN 


Table 4.2 Engine Model Output Parameters 


Variable 

Description 

Symbol 

Part 1 PLA1F_EN 

throttle 1 filtered cockpit power lever angle (deg) 

PLA1FJEN 

PLA2F_EN 

throttle 2 filtered cockpit power lever angle (deg) 

PLA2F_EN 

PLA3F_EN 

throttle 3 filtered cockpit power lever angle (deg) 

PLA3F_EN 

PLA4F.EN 

throttle 4 filtered cockpit power lever angle (deg) 

PLA4F_EN 

Part 2 FXBT_FM 

x-body component of thrust-induced forces (lbf) 

F Xpro P 

FYBT_FM 

y-body component of thrust-induced forces (lbf) 

F Yprop 

FZBT_FM 

z-body component of thrust-induced forces (lbf) 

F Zpnr 

MLBTJFM 

x-body component of thrust-induced moments (ft-lbf) 

My 

*Prop 

MMBT_FM 

y-body component of thrust-induced moments (ft-lbf) 

My 

*Prop 

MNBT_FM 

z-body component of thrust-induced moments (ft-lbf) 

M z 

^ Prop 

FN1_EN 

installed net thrust for engine 1 (lbf) 


FN2_EN 

installed net thrust for engine 2 (lbf) 


FN3_EN 

installed net thrust for engine 3 (lbf) 


FN4_EN 

installed net thrust for engine 4 (lbf) 
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5 Mass Model 


5.1 Mass Model Description 

The Boeing Reference H cycle 1 simulation mass model is discussed in detail in Dawdy et al - 
1994. Mass dynamics are not modeled in the simulation since the mass model does not include 
fuel consumption data. The model contains center of gravity and moment of inertia data for seven 
mass sets. These mass sets are listed in Table 5.1: 


Table 5.1 Mass Model Mass Distributions 


Mass Set 

Description 

Gross Weight 

M01 

operating empty weight 

279,080 (Ibf) 

M05 

M01 + full wing fuel 

646,458 (Ibf) 

M13 

maximum taxi weight at forward eg 

649,914 (Ibf) 

M14 

MZFW + aft body tank full + payload for aft eg 

402,466 (Ibf) 

MIC 

initial cruise condition 

614,864 (Ibf) 

MCR 

mid cruise design point 

501,324 (Ibf) 

MFC 

final cruise condition 

384,862 (Ibf) 


MZFW denotes maximum zero fuel weight. 

5.2 Mass Model Implementation 

The mass model consists of FORTRAN code, which is compiled as a MATLAB FMEX file and 
dynamically linked to the simulation at run time as an S-function. The input parameters to the mass 
model are given in Table 5.2: 


Table 5.2 Mass Model Input Parameters 


Variable 

Description 

Symbol 

AUTOCG FLAG 

eg location flag 


AUTOWT FLAG 

mass set flag 


GW 

gross weight (Ibf) 


GRAV AT 

gravitational acceleration (ft/s'*) 

8 

XCG_PMAC 

longitudinal eg location in percent mean 
aerodynamic chord 


YCG BL MD 

lateral eg wing buttock line location (in) 


ZCG_WL_MD 

vertical eg waterline location (in) 



The input variable AUTOWT_FLAG is used to set the logical variable AUTOWT. If 
AUTOWT_FLAG is equal to 1, AUTOWT is set to true and the mass model buildup equations 
will use data associated with one of the seven mass sets. Gross weight will be set to the gross 
weight for the associated mass set. A zero value for AUTOWTJFLAG causes AUTOWT to be set 
to false. The gross weight input value is used in the table lookups to allow interpolation between 
data points. The input variable AUTOCG_FLAG is used to set the logical variable AUTOCG. If 
AUTOCG_FLAG is set equal to 1, AUTOCG is set to true and the mass model uses the center of 
gravity position associated with the appropriate mass set. An AUTOCG_FLAG value of 0 sets 
AUTOCG to false. Setting AUTOCG to false will cause the mass model to accept the longitudinal 
center of gravity location as a percent of mean aerodynamic chord, lateral center of gravity wing 


24 



buttock line location, and vertical center of gravity waterline location as additional inputs. The 
model uses these values to calculate the center of gravity position; however, it may provide 
inaccurate moment of inertia parameters for a given mass set. The mass model outputs are listed in 


Table 5.3: 


Table 5.3 Mass Model Output Parameters 


Variable 

Definition 

Symbol 

DXCG_MD 

distance from eg to arc along x-body axis (ft) 

X arf 

DYCG_MD 

distance from eg to arc along y-body axis (ft) 

y<2rf 

DZCG.MD 

distance from eg to arc along z-body axis (ft) 

^arf 

WT MD 

gross weight (lbf) 


MASS MD 

mass (slug) 

m 

IXX_MD 

roll moment of inertia (slug-fti) 

h 

IYYJMD 

pitch moment of inertia (slug-ftO 

h 

IZZJMD 

yaw moment of inertia (slug-ft*) 

h 

DCZJVED 

roll-yaw product of inertia (slug-ft') 

Ixz 

XCG_PMAC 

longitudinal eg location in percent mean 
aerodynamic chord 


XCG BS MD 

longitudinal eg body station location (in) 


YCG BL MD 

lateral eg wing buttock line location (in) 


ZCG WL MD 

vertical eg waterline location (in) 
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6 Flight Control System Model 


6.1 Accelerometers 

The simulation performs calculations for axial, lateral, and vertical accelerometers. Accelerometer 
output is in g units. Axial accelerometer output is positive for a forward acceleration, lateral 
accelerometer output is positive for an acceleration to the right, and vertical accelerometer output is 
positive for an upward acceleration! The output for accelerometers positioned at the aircraft center 
of gravity is given by 


V 

_ 1 

u + qw-vr + g sin# 

n y 

v + ru — pw- gcosdsin<j) ► 


8 

— (vv + pv -qu- g cos 6 cos </>) 


( 6 . 1 ) 


For accelerometers positioned at any other location, the accelerations with lever arm effects 
included and the 1 g bias removed are given by 
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( 6 . 2 ) 


where sensor location relative to the center of gravity in the body reference frame, in feet, is given 
by [ x s y s zj. Let x s represent distance forward from the center of gravity, y s represent distance 
to the right of the center of gravity, and z s represent distance below the center of gravity. 
Therefore, sensor location is given by 


r " 

x s 

_ i 
~ 12 

■ ~(BS s - BS C „) 

y s 

K - BL ci ) ■ 

/ \ 

*s. 


{-(WL^W^)] 


where BS, BL, and WL denote locations, in inches, along body station, wing buttock line, 
and waterline respectively. 


6.2 Actuators 


The Reference H control surfaces are hydraulically powered. Actuator dynamics are not modeled 
in the simulation. Inclusion of an actuator model will be part of future model development. 


6.3 Flight Control System 

The simulation models the open-loop, bare airframe and does not include a flight control system. 
Addition of a flight control system model will be part of future model development. 
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7 Atmosphere Model 


7.1 1962 Standard Atmosphere Model 

The simulation uses the 1962 Standard Atmosphere Tables for density and speed of sound. The 
tables contain data in cubic spline format for altitudes ranging from sea level to 240,000 feet. The 

interpolated values for speed of sound, density, and density ratio are calculated as a function of 
altitude. The atmosphere model is written in C and is compiled as a CMEX file. During run time, 
the CMEX file is dynamically linked to the simulation. The atmosphere model is taken from 
LaRCsim (Jackson - 1995). 


7.2 Dry den Turbulence Model 

A simple atmospheric turbulence model, based on the Dryden spectra (Hoblit - 1988), is made 
available in the simulation for the purpose of providing an atmospheric disturbance input to the 
aircraft’s dynamics. Dryden atmospheric turbulence can be simulated in the continuous time 
domain by passing the band-limited white noise signal through the following filters: 


Longitudinal : 


w « 


IZ \ 2Va " 1 
v Ar V K* V-+S 
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Lateral and Vertical : 


w = j-. 



v +s ) 

hva* 




2 


A 


(7.1) 


(7.2) 


Here, V is the vehicle’s total air-relative velocity, c] w are the variance of longitudinal and lateral- 
vertical turbulence, respectively, Lu jW are the spatial scale lengths for longitudinal and lateral- 

vertical turbulence, respectively, and At is the simulation integration time step. The values of <J 2 U W 
and Lu >w used in the simulation are a function of altitude and are, for the case of severe turbulence 
(Johnson et al - 1993), given in Table 7.1. Note that the metric units given in Table 7.1 are 
converted to English units in the implementation of these filters in the simulation. To generate the 
three components of turbulence, three distinct uncorrelated Gaussian white noise sources are used. 
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Table 7. 1 RMS values and Scale Lengths for longitudinal and lateral-vertical 
severe turbulence as a function of altitude. 


Altitude 

(MSL) 

km 

RMS values 

Scale Lengths I 

Longitudinal 
C u , m/s 

Lateral and 
Vertical 
a W3 m/s 

Longitudinal 
U, km 

Lateral and 
Vertical 
Lw, km 

1 

5.70 

4.67 

0.832 

0.624 

2 

5.80 

4.75 

0.902 

0.831 

4 

6.24 

5.13 

1.04 

0.972 

6 

7.16 

5.69 

1.04 

1.01 

8 

7.59 

5.98 

1.04 

0.98 

10 

7.72 

6.00 

1.23 

1.10 

12 

7.89 

5.71 

1.80 

1.54 

14 

6.93 

5.05 

2.82 

2.12 

16 

5.00 

4.31 

3.40 

2.60 

18 

4.07 

3.81 

5.00 

3.34 

20 

3.85 

3.38 

8.64 

4.41 

25 

4.34 

3.34 

12.0 

6.56 

30 

5.60 

3.59 

28.6 

8.88 




8 Simulation Structure 


8.1 Simulation Layout 

Figure 8.1 shows a top level layout of the nonlinear simulation. The simulation is initialized at 
time t = 0 seconds by the initial state vector x 0 . At time t > 0 seconds, the state vector which 
drives the simulation is switched from the initial state vector to the state vector calculated during the 
previous time step. As seen in Figure 8.1, all calculations are performed in nine separate 
subsystems. These subsystems are discussed below. 

Atmospheric disturbances are modeled in th eWind and Turbulence subsystem. This subsystem 
maybe used to specify steady wind velocity components in a local horizontal reference frame and 
atmospheric turbulence velocity components in body axis system, which are used to determine the 
air mass velocity at the aircraft center of gravity. A Dryden atmospheric turbulence model is 
provided in this sub-block which maybe replaced by other atmospheric turbulence models. The 
airmass velocity components are then passed into the Auxiliary Equations subsystem which 
calculates the relative velocity at the aircraft center of gravity. The relative velocity components are 
transferred to the aerodynamic reference center as discussed in Section 3.3. Auxiliary quantities 
needed by the simulation and the aircraft model are computed in this subsystem. The 1962 
Standard Atmosphere model is used in this subsystem. 

The aircraft is modeled by the Mass, Control Inputs, Propulsion, and Aerodynamics subsystems. 
The Mass subsystem executes the Reference H mass model and also calculates the body frame 
gravitational forces acting on the aircraft. The Control Inputs subsystem sets the control surface 
positions needed for the evaluation of the Refere nc e H aerodynamic model and the throttle settings 
for the Reference H engine model. The Reference H engine model is called in the Propulsion 
subsystem. Finally, the aircraft model is completed by the Aerodynamics subsystem which 
executes the Reference H aerodynamic model. 

The External subsystem calculates the total external and inertial forces and moments acting on the 
aircraft at the center of gravity. The propulsive and steady flow aerodynamic forces and moments 
are transferred from the aerodynamic reference center to the center of gravity and added to the 
gravitational force and inertial force components (equation (2.26)) to yield die right-hand-side of 
equation (2.29). The Equations of Motion subsystem solves equation (2.29) for the acceleration 
vector. The six remaining state derivatives are calculated as discussed in Sections 2.6 and 2.7. The 
state derivative vector is then integrated to give the current state vector, which is fed back for use in 
the next simulation time step. 

The Acceleration subsystem calculates accelerometer output at the center of gravity and at a 
specified sensor location. 


8.2 Simulation Implementation 

The SIMULINK implementation of the nonlinear simulation is shown is Figure 8.2. The 
simulati on contains ten subsystems. Each subsystem is discussed below. 

The Initialization subsystem utilizes a switch to control which state vector is used in the simulation. 
At time t = 0 seconds, the initial simulation states are set by the vector InitialStates, which is 
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Figure 8.1 Nonlinear Simulation Layout 
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Figure 8.2 Nonlinear Simulation S1MULINK Implementation - refhsim.m 











defined in the MATLAB workspace. At time t > 0.000001 seconds, the state vector from the 
previous time step is used. Simulation time is also maintained in this subsystem. The 
GustAndWind subsystem uses the Dry den turbulence model to simulate atmospheric disturbances. 
This subsystem calls the functions DrydenVel and dryden_model_coe ff. Steady wind velocities in 
local horizontal coordinate system can be set in this subsystem and atmospheric turbulence velocity 
components in body axes corresponding to severe turbulence levels are calculated. The variable 
ROTATIONJFLAG, which is defined in the MATLAB workspace, is used to include Earth’s 
rotational effects in the equations of motion. If ROTATION_fX,AG is set to 1, Earth rotation is 
included in the computations. If ROTATION_FLAG is set to 0, the computations do not 
include Earth rotation. The variable DRYDENJFLAG is also defined in the MATLAB 
workspace. If DRYDEN_FLAG is set to zero, the Dryden turbulence model is bypassed and 
gust velocities are set to zero. If DRYDEN_FLAG is set equal to 1, then the Dryden turbulence 
model is used to compute translational gust velocities. 

The AuxEqn subsystem calls the function AuxEqn , which computes all necessary quantities at the 
center of gravity and at the aerodynamic reference center. AuxEqn also executes the 1962 Standard 
Atmosphere model function, atmosphere _62. 

The Controls subsystem sets the control positions that are necessary for the Reference H 
aerodynamic and engine models. The control inputs are set by the 22x1 vector control, which is 
defined in the MATLAB workspace. 

The MassModel subsystem executes the Reference H cycle 1 mass model, refhmass, and also 
calls the function GravForce. The function GravForce calculates the body frame components of 
the gravitational force. Gross weight is set by the variable GW, which is defined in the MATLAB 
workspace. The mass model logical variables are set by the variables AUTOWT_FLAG and 
AUTOCG_FLAG, which are also defined in the MATLAB workspace. If AUTOCG_FLAG 
is 0, center of gravity position is set by defining the variables XCG_PMAC, YCG_BL_MD, 
and ZCG_WL_MD in the workspace. 

Engine dynamics are implemented in the EngineModel subsystem. This subsystem also executes 
the engine model refhengine. The engine model logical variable TLIMIT is set by the variable 
TLIMIT_FLAG, which is defined in the MATLAB workspace. 

The AeroModel subsystem executes the Reference H cycle 1 aerodynamic model refhaero. The 
aerodynamic model logical variables are set by the variables INIT_FLAG and RIGIDJFLAG, 
which are defined in the MATLAB workspace. This subsystem calls the function AeroCoejf, 
which computes the steady flow aerodynamic force and moment coefficients at the center of 
gravity. 

The External subsystem calls the function External, which computes the total external and inertial 
forces and moments acting at the center of gravity. The variable ROTATTON_FLAG is used in 
this subsystem. 

The EoM subsystem calls the function Eom .which executes the equations of motion. This 
function calculates the state derivative vector as well as flight path angle, horizontal flight path 
angle, time rate of change of angle of attack at the center of gravity, time rate of change of sideslip 
angle, and time rate of change of total airspeed at the center of gravity. The variable 
ROTATIONJFLAG is also used in this subsystem. 

The last subsystem. Accelerometer, calls the function Accel to compute accelerometer output. 
Sensor location is set by the 3x1 vector sensor, which is defined in the MATLAB workspace. 
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8.3 Simulation Initialization And Setup 


A MATLAB m-file can be used to initialize the aerodynamic model for the nonlinear simulation. 
An example m-file for initialization is given in Figure 8.3. 


% INITIALIZATION 


UAeroTemp = 
1 . ; 

0 ; 

250. ; 

50.; 

- 100 . ; 

0 . 1 ; 

-0.1; 

0 . 1 ; 

50.; 

1 - ; 

- 2 .; 

12 . ; 

- 6 .; 

800. ; 

0.4; 

300. ; 

50.; 

10 . ; 

10 .; 

20 . ; 

10 .; 

10 .; 

30.; 

-10. ; 

30. ; 

30. ; 

20 . ; 

30. ; 

- 20 . ; 

5.; 

-5. ; 

5.; 

5.; 

-15.; 

40. ; 

60. ; 

501324.; 

] ; 


% INIT_FLAG (1 -> yes, 0 -> no) 

% RIGID_FLAG (0 -> elastic, 1 -> rigid) 
% U_ARF (ft/s) 

% V_ARF (ft/s) 

% W_ARF (ft/s) 

% PSA_AR (rad /sec) 

% QSA_AR (rad /sec) 

% RSA_AR (rad /sec) 

% DXCG_MD (ft) 

% DYCG_MD (ft) 

% DZCG_MD (ft) 

% ALPHAA_AR (deg) 

% BETAA_AR (deg) 

% VTOTAL (ft/s) 

% MACH 

% QBAR_AR (psf) 

% HR_ER (ft) 

% DLE1 (deg) 

% DLE2 (deg) 

% DLE3 (deg) 

% DLE4 (deg) 

% DTE1 (deg) 

% DTE2 (deg) 

% DTE 3 (deg) 

% DTE4 (deg) 

% DTE 5 (deg) 

% DTE 6 (deg) 

% DTE7 (deg) 

% DTE 8 (deg) 

% DRUDl (deg) 

% DRUD2 (deg) 

% DRUD3 (deg) 

% DSTAB (deg) 

% DELEV (deg) 

% DVF (deg) 

% DGEAR (deg) 

% GW or WT_MD (lbs) 


YAeroTemp = ref haero (0 ,[], UAeroTemp, 3) ; 


clear UAeroTemp YAeroTemp; 


Figure 8.3 Example Initialization M-file - aero_init.m 


UAeroTemp is an arbitrary input vector used to initialize the aerodynamic model. The only element 
that has any significance is the first element, which is used to set the initialization variable. This 
element must be set equal to 1 . The command line 


YAeroTemp = ref haero ( 0 , [ ] , UAeroTemp , 3 ) ; 


( 8 . 1 ) 
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initializes the aerodynamic model with UAeroTemp by loading all aerodynamic data into memory. 
The initialization process takes approximately 120 seconds. The flight condition of interest is set 
by defining the variables GW, INITJFLAG, RIGID_FLAG, AUTOCG_FLAG, 
AUTOWT_FLAG, TLIMIT_FLAG, ROTATIONJFLAG, DRYDEN_FLAG, 
XCG_PMAC, YCG_PMAC, and ZCG_WL_MD in the MATLAB workspace. The variable 
INIT_FLAG should be set to 0, otherwise the aerodynamic model will be initialized in each time 
step. The vectors InitialStates, control, and sensor must also be defined in the MATLAB 
workspace. Please note that if any of the workspace variables are changed during a 
session, the SIMULINK block diagram must be reinitialized. This reinitialization is 
accomplished by typing the following command line at the MATLAB prompt: 

sizes = refhsim; (8.2) 


An example m-file for setting the simulation is given in Figure 8.4: 

% SET FLIGHT CONDITION FOR YAW DOUBLET DYNAMIC CHECK CASE 


GW = 279080. ; 
INIT_FLAG = 0 . ; 
RIGID_FLAG = 0 . ; 
AUTOWT_F LAG = 1 . ; 
AUTOCG_F LAG = 1 . ; 
TLIMIT_FLAG = 1 . ; 
ROTATION_FLAG = 0. ; 
DRYDEN_F LAG = 0. ; 
XCG_PMAC = 54.602; 
YCG_BL_MD = 0 . ; 
ZCG_WL_MD = 201.5; 



InitialStates = [ 
1036.57435; 

% 

u (ft/s) 

0.; 

% 

v (ft/s) 

27.092043; 

% 

w (ft/s) 

0.; 

% 

p (rad/s) 

0.; 

% 

q (rad/s) 

0.; 

% 

r (rad/s) 

20000. ; 

% 

h (ft) 

0 . *pi/180 . ; 

% 

lat (rad) 

0 . *pi/180 . ; 

% 

Ion (rad) 

-7 . 585 6227 5e-22*pi/ 180 . ; 

% 

phi (rad) 

1 . 497 179 97 5 *pi/ 180 . ; 

% 

theta (rad) 

0 . 0003 19 99 8 *pi/ 180 . ; 
] ; 

% 

psi (rad) 

control = [ 
10 . ; 

% 

DLEl (deg) 

0.; 

% 

DLE2 (deg) 

0. ; 

% 

DLE3 (deg) 

10.; 

% 

DLE4 (deg) 

3. ; 

% 

DTE1 (deg) 

3 - ; 

% 

DTE2 (deg) 

0.; 

% 

DTE 3 (deg) 

0. ; 

% 

DTE 4 (deg) 

0. ; 

% 

DTE 5 (deg) 

0.; 

% 

DTE 6 (deg) 

3.; 

% 

DTE 7 (deg) 

3.; 

% 

DTE 8 (deg) 

0. ; 

% 

DRUD1 (deg) 

0.; 

% 

DRUD2 (deg) 

0.; 

% 

DRUD3 (deg) 

-.426949; 

% 

DSTAB (deg) 
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-.853898; 

0 .; 

0 .; 

68.4916 

68.4916 

68.4916 

68.4916 
] ; 

% Sensor Location 
sensor = [ 

150. ; 

-40. ; 

-36; 

] ; 


% DELEV (deg) 

% DVF (deg) 

% DGEAR (deg) 

% PLA1_CP (deg) 
% PLA2_CP (deg) 
% PLA3_CP (deg) 
% PLA4_CP (deg) 


% 

BS 

sensor 

(in) 

% 

BL 

sensor 

(in) 

% 

WL 

sensor 

(in) 


Figure 8.4 Example Flight Condition Setup M-file - yawdoub_setup.m 


8.4 Simulation Results 

Time history results are given in Appendix B for the flight condition defined in Figure 8.4, with a 
10 degree doublet applied to all three rudder segments. Table 8.1 lists the constants which are hard 
coded into the simulation. In this example simulation, it is assumed that there are no atmospheric 
disturbances and that Earth’s rotation is zero. The simulation was run at 40 Hz using an Euler 
integration algorithm. The simulation took 122.444 seconds to run. These time history results are 
compared with those obtained from independent simulations at the Boeing Commercial Aircraft 
Company and the real-time VMS (Visual Motion Simulator) simulation at NASA Langley Research 
Center. 


Table 8. 1 Constants used in the simulation 


Symbol 

Description 

Numerical Value 

b 

wing span reference length, ft 

129.643 

c 

wing chord reference length, ft 

86.023 

go 

gravitational acceleration at sea level, ft/s 1 

32.228 

R e 

Earth radius, ft 

20,898,908 

£ 

S 

wing reference area, ft 2 

7100.0 

0) E 

Earth rotation, rad/s 

7.2722e-05 
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9 THm Solutions 


9.1 Finding A Trim Solution 

In order to describe the process for finding a trim solution, consider a simple longitudinal level- 
flight trim problem. For the nonlinear simulation, this trim problem can be stated as follows: 

Find 

X = [5, m6 8 e S plal S pla2 8 pla3 8 pla4 u w 6] (9.1) 

such that 

G = r,-/,*, Vr « «]={0} (9.2) 

for a given flight condition. Stabilizer deflection and elevator deflection are denoted by 8 slab and 
8 e respectively. The four throttle positions are defined by 8 plai , 8 pUt2 , 8 pla3 , and S pla4 . These 
throttle positions are equivalent to the cockpit power lever angles discussed in Chapter 4. The 
subscript ‘ des ’ denotes desired value. Independent parameters in equation (9.1) provide the 
necessary degrees of freedom required to satisfy the nonlinear constraints specified in equation 
(9.2), i.e., to perform a longitudinal trim of the aircraft at the desired flight condition. 


9.2 Trim Solution Implementation 

If the elevator is considered to be geared to the stabilizer, then if die stabilizer deflection required 
to trim is known, the corresponding elevator deflection is also known. Therefore, elevator 
deflection required to trim can be removed from in equation (9.1). Also, all four throtties are 
assumed to have equal settings and therefore only one throttle position is required in X . With 
these simplifications, the longitudinal level-flight trim problem can now be stated as follows: 

Find 

X = [(U Splal “ flj (9.3) 

such that 

g=[m-m* s Yv-r.jps V T “ «] ={°} (9- 4 ) 

for a given flight condition. The SIMULINK block diagram used for finding a trim solution is 
shown in Figure 9.1. Engine dynamics and atmospheric disturbances are not modeled. Figure 9.2 
gives an example m-file used to define and solve the longitudinal level-flight trim problem: 
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Figure 9.1 Trim Solution SIMULINK Implementation - refhsim_trim.m 











% SETUP AND SOLVE LONGITUDINAL LEVEL-FLIGHT TRIM PROBLEM 


% SET WORKSPACE VARIABLES 
GW = 501324; 

INIT_FLAG = 0. ; 

RIGID_FLAG = 0 . ; 

AUTOCG_FLAG = 1; 

AUTOWT_FLAG = 1; 

TLIMIT_FLAG = 1 . ; 

ROT AT I ON_F LAG = 0. ; 

XCG_PMAC = 48.0; 

YCG_BL_MD = 0.0; 

ZCG_WL_MD = 226.0; 

ZCG_WL_MD = 190.0; 

% SET SENSOR LOCATION 
sensor = [ 

150. ; % BS (in) 

-30.; % BL (in) 

-36. ; % WL (in) 

] ; 

% REINITIALIZE BLOCK DIAGRAM 
sizes = refhsim_trim; 

% SET CONTROL INPUTS (minus DSTAB, DELEV, and PLA(i)_CP) 


U = [ 

10. ; 

% 

1 

- DLEl (deg) 

0. ; 

% 

2 

- DLE2 (deg) 

0. ; 

% 

3 

- DLE3 (deg) 

10. ; 

% 

4 

- DLE4 (deg) 

3. ; 

% 

5 

- DTE1 (deg) 

3. ; 

% 

6 

- DTE 2 (deg) 

0. ; 

% 

7 

- DTE3 (deg) 

0. ; 

% 

8 

- DTE4 (deg) 

0. ; 

% 

9 

- DTE 5 (deg) 

0. ; 

% 

10 

- DTE 6 (deg) 

3. ; 

% 

11 

- DTE7 (deg) 

3. ; 

% 

12 

- DTE 8 (deg) 

0. ; 

% 

13 

- DRUD1 (deg) 

0. ; 

% 

14 

- DRUD2 (deg) 

0. ; 

% 

15 

- DRUD3 (deg) 

0 • ; 

% 

16 

- DVF (deg) 

0. ; 
]? 

% 

17 

- DGEAR (deg) 

% SET STATES 
x = [ 

(minus u, 

w, 

and theta) 

0. ; 

% 

1 - 

v (ft/s) 

0. ; 

% 

2 - 

p (rad/s) 

0. ; 

% 

3 - 

g (rad/s) 

0. ; 

% 

4 - 

r (rad/s) 

30000. ; 

% 

5 - 

h (ft) 

0 . *pi/18Q . ? 

% 

6 - 

lat (rad) 

0 . *pi/180 . ; 

% 

7 - 

Ion (rad) 

0 . *pi/180 . ; 

% 

8 - 

phi (rad) 

-140 . *pi/180 . ; 
] ; 

% 

9 - 

psi (rad) 

% SET DESIRED 

CONSTRAINT 

VALUES 

g des = [ 




0.8; 

% 

Mach 

0 . 0*pi/180 ; 

% 

gamma 

0.0; 

% 

VT 

_dot 
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0 . 0*pi/180 ; 
0 . 0*pi/180 ; 
] ; 


% alpha_dot 
% q_dot 


% SET OPTIONS VECTOR FOR CONSTR FUNCTION 

% 12 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 

options = [ 1 le-4 le-4 le-6 00000 0 0 0 5 1000 0 le-6 1.0 1.0]'; 

% SET INITIAL GUESS, LOWER BOUNDS, AND UPPER BOUNDS for X 


% 

DSTAB 

PLA CP 

u 

w 

Theta 

Xstart = [ 

3. 

50 

900 

65. 

2 . 0*pi/180 ] ' ; 

Xlb = [ 

-15 

0 

0 

-500 

-15*pi/180 ] ' ; 

Xub = [ 

15 

100 

2600 

500 

35*pi/180 ]'; 


% FIND TRIM SOLUTION 

Xfinal = constr( ' longtrim_fun' ,Xstart, options, Xlb.Xub, ' ' ,u,x,g_des) 

% COMPUTE OBJECTIVE FUNCTION, CONSTRAINT VECTOR, AND OUTPUT VECTOR 
% AT TRIMMED CONDITIONS 
[ f , G, Y] =longtrim_fun (Xfinal , u, x, g_des ) 

% DISPLAY RESULTS 

fprintf (’ Longitudinal , Level-Flight Trim Solution\n\n'); 
fprintf ( 1 DSTAB = %10.5f DELEV = %10.5f PLA_CP=%10 . 5f \n ' , ... 

Xfinal (1) , 2*Xfinal (1) , Xfinal (2) ) ; 

f print f ( 'u = %10 . 5f w = % 1 0 . 5 f theta = %10.5f\n', ... 

Xfinal (3) , Xfinal (4) , Xfinal (5) *180 /pi) ; 
fprintf (’M = %10.5f gamma = %10.5f VT_dot = %10.5f\n', ... 

Y (3 ) , Y (7 ) *180 /pi , Y (23 ) ) ; 

fprintf (' alpha_dot = %10.5f qdot = %10 . 5f \n‘ , Y(21) ,Y(15) ) ; 

Figure 9.2 Example M-File For Computing A Longitudinal Trim Solution - longtrim.m 


Please see the Optimization Toolbox User's Guide for a detailed explanation of the vector 
options. The command line 


sizes = refhsim_trim; (9.5) 

reinitializes the SIMULINK block diagram as discussed in Section 8.2. A constrained nonlinear 
optimization is performed by the MATLAB function constr on the function longtrim_fun. The 
function file longtrim_fun.m is given in Figure 9.3. The outputs of the trim process are the 
elements of the trim input vector X , the block diagram output vector Y , the constraint vector G , 
and the objective function f. 

function [ f , G, Y] = longtrim_fun(X,u,x,g_des) 

% FIND TRIM SOLUTION 

% SET INITIAL CONDITIONS 
xO = [ ] ; 

uO = [ 
u(l) ; 
u(2) ; 
u ( 3 ) ; 
u (4) ; 
u(5) ; 
u(6) ; 

U(7) ; 
u(8) ; 
u(9) ; 
u(10) ; 


% 1 - DLE1 (deg) 
% 2 - DLE2 (deg) 
% 3 - DLE3 (deg) 
% 4 - DLE4 (deg) 
% 5 - DTE1 (deg) 
% 6 - DTE2 (deg) 
% 7 - DTE3 (deg) 
% 8 - DTE4 (deg) 
% 9 - DTE 5 (deg) 
% 10 - DTE6 (deg) 
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u(ll) ; 
u(12) ; 
u(13 ) ; 
u(14); 
u(15); 
X(l) ; 
2*X(1) ; 
u ( 1 6 ) ; 
u ( 17 ) ; 
X(2) ; 

% 

% 

% 

X(3) ; 
x(l) ; 
X(4) ; 
x (2 ) ; 
x(3) ; 
x(4) ; 
x ( 5 ) ; 
x ( 6 ) ; 
x(7) ; 
x(8) ; 
X(5) ; 
x(9) ; 

] ; 


% 11 - DTE7 (deg) 

% 12 - DTE8 (deg) 

% 13 - DRUD1 (deg) 

% 14 - DRUD2 (deg) 

% 15 - DRUD3 (deg) 

% 16 - DSTAB (deg) 

% 17 - DELEV (deg) -> DELEV = 2* DSTAB 
% 18 - DVF (deg) 

% 19 - DGEAR (deg) 

% 20 - PLA1_CP (deg) 

- PLA2_CP (deg) 

- PLA3_CP (deg) 

- PLA4_CP (deg) 

% 21 - u (ft/s) 

% 22 - v (ft/s) 

% 23 - w (ft/s) 

% 24 - p (rad/s) 

% 25 - g (rad/s) 

% 26 - r (rad/s) 

% 27 - h (ft) 

% 28 - lat (rad) 

% 29 - Ion (rad) 

% 30 - phi (rad) 

% 31 - theta (rad) 

% 32 - psi (rad) 


% CALCULATE OUTPUTS 
Y = refhsim_trim(0,x0,u0, 3) ; 


% CALCULATE CONSTRAINTS AND 
G ( 1 ) = Y ( 3 ) - g_des ( 1 ) ; 

G (2 ) = Y(7) - g_des(2); 

G (3 ) = Y(23 ) - g_des ( 3 ) ; 

G(4) = Y (21 ) - g_des ( 4 ) ; 

G(5) = Y(15) - g_des ( 5 ) ; 

f = 0.0; 


OBJECTIVE FUNCTION 
% Mach Number constraint 
% Flight Path angle constraint 
% VT_dot constraint 
% alpha_dot constraint 
% q_dot constraint 


Figure 9.3 Longitudinal Trim Function longtrim_fun - longtrim_fun.m 
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10 Linear Models 


10.1 Introduction 

This chapter describes the methodology by which linear models are generated from the nonlinear 
simulation. Linear model of the aircraft dynamics can be obtained from the nonlinear simulation in 
MATLAB/SIMULENK environment by using the MATLAB’s linmod command. However, since 
the nonlinear simulation uses the body axes equations of motion, the linear models generated are 
body axes linear models, which are typically not convenient for controls engineering. Therefore, 
the body axes state-space linear models are transformed to a set of preferred states and outputs. 
This chapter describes that transformation process and presents an example case. 


10.2 States, Inputs, and Outputs 

Let the nonlinear simulation state vector, input vector, and output vector be given by x, u, and 
y respectively. Also, let X define a new state vector with the preferred states, u redefines the 
input vector, and y define a new output vector with the preferred outputs. The state vectors, 
input vectors, and output vectors previously mentioned are defined in Table 10. 1. For the engine 
dynamics, t ll represents the lead-lag filter state and ( 0 B represents the lag filter state. y h denotes 
the horizontal flight path angle, y v denotes the vertical flight path angle, n y denotes the lateral 
acceleration at the center of gravity, n z denotes the vertical acceleration at the center of gravity, n y ^ 
denotes the lateral acceleration at the pilot station, and n z ^ denotes the vertical acceleration at the 
pilot station. 


10.3 State Vector and Output Vector Transformations 

Looking at the two state vectors, the states V T , a, and /? replace the states u, v, and w. 

Recall that total airspeed, angle of attack, and sideslip angle are functions of the body frame 
components of inertial translational velocity; therefore, the two state vectors are related as follows: 

x = T x x, (10.1) 

where T x is some undefined transformation matrix. Similarly, the output vectors are related by 

y = T y y , (10.2) 

where T y is another undefined transformation matrix. To define T x and T y , recall that total airspeed is 
given by 


V T = V « 2 + v 2 +w 2 . 

Expanding V T about V Tq using Taylor series gives 


(10.3) 
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AV t = V Tu Au + ^V Tuu (Au) 2 +... + V Tv Av + ^V Tvv (Av ) 2 + ... + V Tw Aw + ^V Tww (Aw) 2 +... (10.4) 


Table 10.1 State, Input, and Output Vectors. 


Element 

X 

u 

y 

X 

u 

y 

1 

4 

DLE1 

U 

v T 

DLE1 

Vt 

2 

X LL 4 

DLE2 

V 

a 

DLE2 

a 

3 

a> B3 

DLE3 

w 

q 

DLE3 

q 

4 

t LL 3 

DLE4 

p 

e 

DLE4 

e 

5 

C0 B2 

DTE1 

q 

h 

DTE1 

h 

6 

T LL 2 

DTE2 

r 

P 

DTE2 

P 

7 

®B1 

DTE3 

h 

r 

DTE3 

r 

8 

Z LL 1 

DTE4 

0 

0 

DTE4 

0 

9 

u 

DTE5 

e 

P 

DTE5 

P 

10 

V 

DTE6 

¥ 

¥ 

DTE6 

¥ 

11 

w 

DTE7 

z LLl 

A 

DTE7 

t ll\ 

12 

p 

DTE8 

a B\ 

T 

DTE8 

®B 1 

13 

q 

DRUD1 

T LL 2 

t LL 1 

DRUD1 

Z LL 2 

14 

r 

DRUD2 

(0 B2 

®BI 

DRUD2 

(0 B 2 

15 

h 

DRUD3 

z U.‘i 

r LL 2 

DRUD3 

r LL 3 

16 

A 

DSTAB 

®B 3 

®B2 

DSTAB 

co B3 

17 

T 

DELEV 

t lla 

T ZX3 

DELEV 

T LL 4 

18 

0 

DVF 

a B 4 

®B3 

DVF 

0) B4 

19 

0 

DGEAR 

X/, 

r LL 4 

DGEAR 

7h 

20 

¥ 

PLAl.CP 

Xv 

®B4 

PLA1_CP 

Xv 

21 


PLA2.CP 

n y 

PLA2_CP 

n y 

22 


PLA3_CP 

n z 


PLA3_CP 

n z 

23 


PLA4_CP 

% 


PLA4.CP 

% 

24 



\ 




25 



h 



h 

26 



% 



% 

27 






% 

28 



F z 

Z *f 



% 

29 






M x, 

30 



M r< 



M r, 

31 



M z, 



M Zsf 


Keeping only the linear terms from equation (10.4) yields 


AVr = V t Au + V t Av + V t Aw, 

1 l u Ay A W 


(10.5) 


where 


dVj 

du 


x=x 0 


(10.6,a) 
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(10.6, b) 


V T =& 


V T = 


dv 

dv T 


ix=x 0 


dw 


(10.6,c) 


Equation (10.5) relates the linear variation in total airspeed to the linear variations in inertial 
translational velocity body frame components. If the procedure just described is applied to angle of 
attack and sideslip angle, the linear variation in angle of attack about a 0 is given as 


A a = a u Au + a w Aw 
and the linear variation in sideslip angle about j3 0 is given by 

A/3 = P U A u + /3 v Av + ft, Aw. 


(10.7) 


( 10 . 8 ) 


Now the transformation matrices can be defined. The transformation matrix T x is given by 


where 


Ts = 


T — 

®12x8 

TSi2xl2 

A x ” 

OO 

X 

OO 

£ 

1 

®8xl2 . 


V T 

Vj 

A v 

V T 

0 

0 

0 

0 

0 

0 

0 

0 

0 

ct u 

0 

a w 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

ft 

ft 

Pw 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 


(10.9) 


( 10 . 10 ) 


and 
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Te = 


( 10 . 11 ) 


0 0 0 0 0 0 0 f 
0 0 0 0 0 0 1 0 
0 0 0 0 0 1 0 0 
0 0 0 0 1 0 0 0 
0 0 0 1 0 0 0 0 
0 0 1 0 0 0 0 0 
0 1 0 0 0 0 0 0 

1 0 0 0 0 0 0 0 


The elements of Ts which are functions of inertial body frame acceleration are given by 


u 

v T =— 
v T 


V T =■ 


V T = — 
* V T 

—w 


a u ~ — t 


u 


1 + 1 - 


2^\ 


a w =- 


Pu=- 


v j 

“MV 


vl 


1- 


( ^ 
I v * 




Pv=- 


-1 


1- 


-T 

\Vt) 


K V T Vfj 


Av=- 


— VW 



V Tj 


\ 2 


(10. 12, a) 
(10.12,b) 
(10.12,c) 
(10.12,d) 

(10.12, e) 

(10.12,f) 

(10.12,g) 

(10.12, h) 
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The transformation matrix T y is given by 


where 


To = 



T — 

To 

10x10 

® 10x21 


h ~ 

L ®21xl0 

^21x21 

V T 

V T 

*v 

V? 

‘ H’ 

0 

0 0 

a u 

0 


0 

0 0 

0 

0 

0 

0 

1 0 

0 

0 

0 

0 

0 0 

0 

0 

0 

0 

0 0 

0 

0 

0 

1 

0 0 

0 

0 

0 

0 

0 1 

0 

0 

0 

0 

0 0 

Pu 

Pv 

p w 

0 

0 0 

0 

0 

0 

0 

0 0 


0 


(10.13) 


(10.14) 


10.4 State-Space Form 

Let the linear model generated from the nonlinear simulation be given in state-space form by 


x = Ax + Bu 
y = Cx + Du . 

Let the transformed linear model with the desired states and outputs be given by 

x = Ax + Bu 
y = Cx + Du. 

This section describes how equation (10.15) is transformed to give equation (10.16). 
Leftmultiplying both sides of equation (10.1) by T x 1 gives 

T“*x = x , 

and taking the derivative with respect to time of equation (10.17) yields 

T x _, x = x. 


(10.15,a) 

(10.15,b) 

(10. 16, a) 
(10.16, b) 

(10.17) 

(10.18) 


Substituting equations (10.17) and (10.18) into equation (10.15,a) and rearranging yields 
equation (10. 16, a), where 

A = T X AT X 1 (10. 19, a) 
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and 


B = T X B. (10. 19,b) 

Following the process outlined above, equation (10.16,b) is obtained, where 

C = TyCT" 1 (10.20,a) 

and 

D = T y D. (10.20,b) 


10.5 Units 

The simulation states and outputs have units of radians for all angular quantities. Therefore, the 
linear models also have units of radians. For convenience, let the linear models be expressed in 
units of degrees. To convert the linear models generated from the simulation from units of radians 
to units of degrees, let equation (10.16) be rewritten as 

x r = A r x r + h r u d (10.21, a) 

y r =C r x r + D r u^, (10.21,b) 

where the symbols * r ’ and * d ' denote radians and degrees respectively. The state vector in units 
of degrees is related to the state vector in units of radians by 


where 


x d -S x x r< 


( 10 . 22 ) 


SSl2xl2 

®12x8 

. ®8xl2 

1 8x8 . 


(10.23) 


where 



Ss = 


0 

180 


n 

0 0 

0 0 

0 0 

0 0 

0 0 

0 0 

0 0 

0 0 

0 0 

0 0 


0 

0 

180 

tc 

0 

0 

0 

0 

0 

0 

0 

0 

0 


0 0 
0 0 


7T 

0 


0 0 
0 


0 

0 


0 0 0 
180 


0 0 
1 


K 

0 0 


0 0 0 


0 0 0 


0 0 0 


0 0 0 


0 0 0 


0 

0 


0 

180 


0 

0 


0 

0 

0 

180 

K 

0 

0 

0 

0 

0 

-1 


Leftmultiplying both sides of equation (10.22) by S x gives 

Sx\/ = x r . 

Taking the time derivative of equation (10.25) yields 

S; 1 x rf =x r . 


0 

0 


0 

0 

0 

0 

180 

it 

0 

0 

0 

0 


0 

0 


0 0 0 0 


0 

0 

0 

0 

0 

180 

It 

0 

0 

0 


0 

0 


0 

0 

0 

0 

0 

0 

180 

7t 

0 

0 


0 

0 

0 

0 

0 

0 

0 

180 

K 

0 


0 

0 


0 0 


0 

0 

0 

0 

0 

0 

0 

0 

180 
7T . 


(10.24) 


(10.25) 


(10.26) 


Substituting equations (10.25) and (10.26) into equation (10.21,a) and rearranging gives the 
following equation: 


x d =A d x d + B d u d , (10.27) 

where 

A d =S x A r S~ 1 (10.28, a) 

and 

B d =S x B r . (10.28, b) 

Similarly, the following expression relates the output vector in units of degrees to the output 
vector in units of radians: 


Id =s y y r , 


(10.29) 
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where S y is given by 


where 


and 


Ss = 


1 

0 

0 

0 

0 

0 

0 

0 

0 

0 


Sv = 


0 

180 

K 

0 

0 

0 

0 

0 

0 

0 

0 


SSiOxlO 

®10x8 

®10x2 

®10xll 

®8xl0 

^8x8 

®8x2 

®8xll 

®2xl0 

®2x8 

Se 2 x2 

®2xll 

_0llxl0 

®llx8 

®llx2 

lxl I _ 

0 0 

0 

0 

0 0 

0 0 

0 

0 

0 0 

80 

— 0 

0 

0 

0 0 

C 

oc 

0 — 

0 

0 

0 0 

n 
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1 

0 

180 

0 0 

0 0 

0 

K 

0 0 


0 

0 

0 

0 


0 

0 

0 

0 


Se = 


0 

0 

0 

0 

180 

7Z 

0 


0 

0 

0 

0 


180 

n 

0 

0 
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0 
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n 


0 
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n 

0 
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0 

0 

0 
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n 
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0 

0 

0 

0 

0 

0 

0 

0 

0 

180 
n J 


(10.30) 


(10.31) 


(10.32) 


Following the process outlined above, the linear model for the output vector in degrees is given by 

y d = C d x d +'D d u d , (10.33) 

where 


C d — SyCjS 


-1 


(10.34,a) 


and 
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\y d ~s y D r . 


(10.34,b) 


10.6 Linear Model Implementation 

MATLAB uses the function linmod to extract a linear state-space model of a system about an 
operating point, linmod perturbs each state and each input about the operating point to determine 
changes in the state derivatives and the outputs. The states are perturbed by 

x(i) + Ax(i), (10.35) 

where A x(i) is the perturbation associated with the f h state. Similarly, the inputs are perturbed by 

u(j) + Au(j) , (10.36) 

where A u(j) is the perturbation for the f h input. This method only considers the data on “one 
side” of the operating point and may produce incorrect linear models. For example, this 
perturbation method generates incorrect results for the effects of lateral control surfaces on 
longitudinal state derivatives and outputs. Because of this deficiency in linmod, the linmod source 
file was modified to include positive and negative perturbations about an operating point. The 
states are now perturbed by 


x(i)± Ax(i), (10.37) 

and the inputs are perturbed by 

u(J) ± AmO') . (10.38) 

The function linearmodel generates a state-space linear model using the perturbation technique 
given by equation ( 10.37) and equation ( 1 0.38). The input to linearmodel is the operating point 
which is defined by the state vector, X, and the input vector, u . linearmodel outputs the state- 
space linear model given by equation (10.15). 

The SIMULINK block diagram used for generating linear models is shown in Figure 10.1. The 
Dryden Turbulence Model is not included. The lead-lag filter and lag filter in the engine dynamics 
are implemented in state-space form so that the engine states may be easily extracted. Please 
note that any changes to the block diagram may alter the order of the states and 
cause the generation of incorrect linear models. As discussed in Section 8.2, the 
aerodynamic model must be initialized before generating linear models and if the workspace 
variables are changed during a session the block diagram must be reinitialized, as discussed in 

Section 8.2. Figure 10.2 gives an example m-file used to set up the linear model problem. The m- 
file refhlinearmodel.m generates a linear model with the desired states and outputs as given by 
equations (10.27) and (10.33). 


49 



Figure 10.1 Linear Model SIMULINK Implementation - refhsim lin.m 






% GENERATES LINEAR MODEL FOR fc2deg -> Supersonic Cruise Mid 


GW = 0 . 501324e6 ; 
INIT_FLAG = 0. ; 
RIGID_FLAG = 0 . ; 
AUTOWT_FLAG = 1 . ; 
AUTOCG_FLAG = 1 . ; 
TLIMIT_FLAG = 1 . ; 
ROTATION_FLAG = 0 . ; 
XCG_PMAC = 54.8146; 
YCG_BL_MD = 0 . ; 
ZCG_WL_MD = 203.15; 


% SENSOR LOCATION 


sensor = [ 
150. ; 

-30.; 

-36.; 

] ; 


% BS (in) 
% BL (in) 
% WL (in) 


% DEFINE INITIAL STATE VECTOR 
xO = ( 

% ENGINE STATES 


40. ; 

% 

wb4 

(deg/s) 

100. ; 

% 

tll4 

(s) 

40. ; 

% 

wb3 

(deg/s) 

100. ; 

% 

tll3 

(s) 

40. ; 

% 

wb2 

(deg/s) 

100. ; 

% 

tll2 

(s) 

40. ; 

% 

wbl 

(deg/s) 

100. ; 

% 

till 

(s) 

% SIMULATION STATES 



2318.8; 

% 

u (ft/s) 

0. ; 

% 

V (ft/s) 

145.96; 

% 

w (ft/s) 

0. ; 

% 

p (rad/s) 

0. ; 

% 

q (rad/s) 

0. ; 

% 

r (rad/s) 

57794. ; 

% 

h (ft) 

0 . *pi/180 . ; 

% 

lat 

(rad) 

0 . *pi/180 . ; 

% 

Ion 

(rad) 

0. *pi/180. ; 

% 

phi 

(rad) 

3 . 6018*pi/180 . ; 

% 

theta (rad) 

-90 . *pi/180 . ; 
1 ; 

% 

psi 

(rad) 

% DEFINE INITIAL INPUT VECTOR 


uO = [ 

O.OOOOOOOOOOE+OO 

% 

DLEl 

(deg) 

O.OOOOOOOOOOE+OO 

% 

DLE2 

(deg) 

O.OOOOOOOOOOE+OO 

% 

DLE3 

(deg) 

O.OOOOOOOOOOE+OO 

% 

DLE4 

(deg) 

O.OOOOOOOOOOE+OO 

% 

DTE1 

(deg) 

O.OOOOOOOOOOE+OO 

% 

DTE 2 

(deg) 

O.OOOOOOOOOOE+OO 

% 

DTE 3 

(deg) 

O.OOOOOOOOOOE+OO 

% 

DTE4 

(deg) 

O.OOOOOOOOOOE+OO 

% 

DTE 5 

(deg) 

O.OOOOOOOOOOE+OO 

% 

DTE 6 

(deg) 

O.OOOOOOOOOOE+OO 

% 

DTE 7 

(deg) 

O.OOOOOOOOOOE+OO 

% 

DTE 8 

(deg) 

O.OOOOOOOOOOE+OO 

% 

DRUD1 (deg) 

O.OOOOOOOOOOE+OO 

% 

DRUD2 (deg) 

O.OOOOOOOOOOE+OO 

% 

DRUD3 (deg) 
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0. 1892955855 E+01; 
0.3785911710E+01; 
O.OOOOOOOOOOE+OO; 
O.OOOOOOOOOOE+OO; 
0 . 1000000000E+03 ; 
0.1000000000E+03; 
0.1000000000E+03; 
0.1000000000E+03; 
] ; 


% DSTAB (deg) 

% DELEV (deg) 

% DVF (deg) 

% DGEAR (deg) 

% PLA1_CP (deg) 
% PLA2_CP (deg) 
% PLA3_CP (deg) 
% PLA4_CP (deg) 


% STATE-SPACE REPRESENTATION FOR LEAD-LAG FILTER 
ATau = [-4 . ] ; 

BTau = [1. ] ; 

CTau = [-.33333; 1. ] ; 

DTau = [ . 33333; 0] ; 


% STATE-SPACE REPRESENTATION FOR LAG FILTER 
AOmega = [-2.5]; 

BOmega = [ 1 . ] ; 

COmega = [2.5; 1 . ] ; 

DOmega = [ 0 . ; 0 . ] ; 


Figure 10.2 Example M-file To Setup The Linear Model Problem - fc2deg_setup.m 


10.7 Linear Model Results 

Time history results are shown in Figure 10.3 for the linear model generated at the flight 
condition defined in Figure 10.2, with a one degree doublet applied to the horizontal stabilizer 
and a corresponding two degree doublet applied to the elevator. Atmospheric disturbances are 
neglected. Results generated by the MATLAB/SIMULINK nonlinear simulation are also 
presented for comparison. Both sets of time histories were generated at 40 Hz. 
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Figure 10.3 Linear Model And Nonlinear Simulation Dynamic Response To A Stabilizer Doublet 
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Figure 10.3 (Continued) 


54 






aw Rate, BA (deg/s) 













— Nonlinear Simulation 
Linear Model 












Elevator (deg) Stabilizer (deg) 


4 


3 \- 


— Nonlinear Simulation 
Linear Model 


2 1 — 


10 15 

Time (s) 


20 



Figure 10.3 (Concluded) 
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11 Concluding Remarks 


The mathematical model and associated code used to simulate the Boeing Reference H 
configuration high speed civil transport aitcraft are described in this report. The simulation is 
implemented using MATLAB/SIMULINK software, a commercially available control system 
engineering tool, to support advanced control law research. The simulation software includes the 
capability to calculate trim solutions and to generate linear models at specified flight conditions. 

The simulation math model includes the nonlinear, six degree-of-freedom equations which govern 
the motion of a rigid aircraft in atmospheric flight. The Earth atmosphere is modeled with the 1962 
Standard Atmosphere and a Dryden turbulence model. The Reference H configuration is modeled 
using the Boeing Cycle 1 simulation data base, which consists of geometry information, tabular 
aerodynamic model, an engine model, and a mass model. Models of actuator dynamics, landing 
gear, and flight control system are not implemented in the simulation. 

To validate the MATLAB/SIMULINK simulation software, a comparison of simulation-generated 
dynamic responses with those from independent simulations at Boeing Commercial Aircraft 
Company and NASA Langley Research Center’s VMS were made. A good match of the 
responses was achieved for the dynamic check cases provided. The trim solutions obtained with 
the developed software compared closely to the trim check cases provided by Boeing. The linear 
models obtained at representative flight conditions compared closely to those generated using the 
NASA Langley Research Center’s real time simulation used for the Visual Motion Simulation. 
Furthermore, time histories generated using the linear model compared closely to those generated 
by the nonlinear simulation. 

The source code is available from the authors. 
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Symbols and Abbreviations 

Acronyms 


arc aerodynamic reference center 

BA body axis reference frame 

BL location along wing buttock line, in 

BS location along body station, in 

eg center of gravity 

SA stability axis reference frame 

WL location along waterline, in 

Symbol Description 


A,B,C,D 

a 

b 

C L 

Cl 

C n i 

Cy 

C 

{F} 

{F}, 

F 

F u (s ) 
F w (s) 

G 

8 
8 o 

h 

h 

K 

h 

K 

K 

K 

W] A 

[m], 

M 


state-space matrices in the linear model description 

speed of sound, ft/s 

wing span reference length, ft 

drag coefficient 

lift coefficient 

rolling mo ment coefficient 

pitching moment coefficient 

yawing moment coefficient 

sideforce coefficient 

wing chord reference length, ft 

generalized force vector 

generalized inertial force vector 

force, lbf 

longitudinal turbulence frequency response function 

lateral-vertical turbulence frequency response function 

trim constraint vector 

gravitational acceleration, ft/s 2 

gravitational acceleration at sea level, ft/s 2 

altitude above sea level, ft 

roll moment of inertia, slug-ft 2 

roll-yaw product of inertia, slug-ft 2 

pitch moment of inertia, slug-ft 2 

yaw moment of inertia, slug-ft 2 

longitudinal turbulence spatial scale length, ft 

lateral-vertical turbulence spatial scale length, ft 

apparent mass matrix 

inertial mass matrix 

Mach number or 

moment, ft-lbf 

mass, slugs 

axial accelerometer output, g units 
lateral accelerometer output, g units 
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"z 

p 


p 

Q 

q 

q 

R 


Re 

r 

S 

S, 

S s 

S x 

s, 

[ t bv ] 
[^bs] 

T, 

T„ 

T s 

T x 

T , 

t 

u 

u 

{v} 

V D 

V E 

y > ec 
V N 

V T 


X 

x 


y 

y arf 

y* 


Z s 


vertical accelerometer output, g units 

x-body frame component of rotational velocity of the body reference frame relative 

to the local reference frame, rad/s 

x-body frame component of rotational velocity, rad/s 

y-body frame component of rotational velocity of the body reference frame relative 

to the local reference frame, rad/s 

y-body frame component of rotational velocity, rad/s 

dynamic pressure, lbf/ft 2 

distance from center of Earth to aircraft mass center, ft or 

z-body frame component of rotational velocity of the body reference frame relative 

to the local reference frame, rad/s 

Earth radius, ft 

z-body frame component of rotational velocity, rad/s 
wing reference area, ft 2 

sub-matrix of the output units transformation matrix 
sub-matrix of the state and output units transformation matrix 
state units transformation matrix 
output units transformation matrix 

transformation matrix from vehicle-carried local horizontal reference frame to body 
reference frame 

transformation matrix from stability reference frame to body reference frame 

sub-matrix of the state vector transformation matrix 

sub-matrix of the output vector transformation matrix 

sub-matrix of the state vector transformation matrix 

state vector transformation matrix 

output vector transformation matrix 

time, s 
control vector 

x-body frame component of inertial translational velocity, ft/s 
generalized velocity vector 
downward velocity component, ft/s 
eastward velocity component, ft/s 

eastward velocity component corrected for Earth rotation, ft/s 
northward velocity component, ft/s 
total airspeed, ft/s 

y-body frame component of inertial translational velocity, ft/s 
z-body frame component of inertial translational velocity, ft/s 
trim input vector 
state vector 

distance from eg to arc along x-body axis, + forward, ft 
distance from eg to sensor along x-body axis, + forward, ft 
output vector 

distance from eg to arc along y-body axis, + right, ft 
distance from eg to sensor along y-body axis, + right, ft 
distance from eg to arc along z-body axis, + down, ft 
distance from eg to sensor along z-body axis, + down, ft 
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a 

p 

At 

$pla 

^stab 

<t> 

Yh 

7v 

A 

e 

p 

<*u 

<j w 

T 

*LL 

[Q] 

(0 B 

(0 E 

¥ 


angle of attack, deg or rad 
sideslip angle, deg or rad 
integration time step, s 
elevator deflection, deg 
power lever angle, deg 
stabilizer deflection, deg 
Euler roll angle, rad 

horizontal flight path angle measured from north, rad 

vertical flight path angle, rad 

latitude, rad 

Euler pitch angle, rad 

air density, slug/ft 3 

variance of longitudinal turbulence, ft/s 

variance of lateral-vertical turbulence, ft/s 

longitude, rad 

engine lead-lag filter state, s 

body frame inertial rotation matrix 

engine lag filter state, deg/s 

Earth rotation, rad/s 

Euler yaw angle, rad 


Subscripts and Superscripts 


0 

Aero 

a 

airmass 

arf 

b 

E 

d 

Grav 

Prop 

ps 


s 

Sf 

usf 

u 

w 


initial 

due to air loads 

aerodynamic, computed at center of gravity with respect to air mass 

due to winds, gusts, and turbulence 

computed at aerodynamic reference center 

body reference frame 

Earth 

degrees 

due to gravity 

thrust-induced 

pilot station 

radians 

sensor 

steady flow 

unsteady flow 

longitudinal atmospheric turbulence quantities 
lateral- vertical atmospheric turbulence quantities 
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Simulation Files 


Function Files 


File 

Description 

Accel.m 

AeroCoeff.m 

AuxEqn.m 

DrydenVel.m 

dryden_model_coeff.m 

Eom.m 

Extemal.m 

GravForce.m 

linearmodel.m 

longtrim_fun.m 

performs accelerometer output calculations 

transfers steady flow aerodynamic force and moment coefficients 

from the aerodynamic reference center to the center of gravity 

performs auxiliary calculations which drive the aircraft 

model and the equations of motion 

calculates airspeed for use in the Dryden turbulence model 

calculates turbulence frequency response function coefficients 

calculates the state derivative vector and auxiliary quantities 

calculates total external forces and moments acting at the 

center of gravity 

calculates gravitational force components 
generates linear model in state-space form 
trim optimization function 

S-function MEX Files 


atmosphere_62 .mex4 
refhaero.mex4 
refhengine.mex4 
refhmass.mex4 

1962 Standard Atmosphere model 
Reference H cycle 1 aerodynamic model 
Reference H cycle 1 engine model 
Reference H cycle 1 mass model 

SIMULINK Files 


refhsim.m 

refhsim_lin.m 

refhsim_trim.m 

nonlinear simulation block diagram 
linear model block diagram 
trim solution block diagram 

Auxiliary M-Files 


aero_init.m 

fc2deg_setup.m 

longtrim.m 

plot_var.m 

refh_linearmodeI.m 

yawdoub_setup.m 

initializes aerodynamic model 

sets up flight condition and simulation to compute linear model at the 

desired flight condition. 

sets up and finds longitudinal trim solution 

generates plotting variables from simulation output 

first computes the linear model in body axis states and then 

transforms linear model to desired states and output 

sets initial flight condition and simulation variables for rudder 

doublet dynamic response 
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Appendix A - Apparent Mass Matrix Elements 

Each element of the 6x6 apparent mass matrix is defined below. Given the following two 
intermediate terms: 
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The apparent mass matrix elements are given by 
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(A. 12) 
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(A. 17) 
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The elements given by equations (A.3) to (A. 1 1) have units of slugs. Elements (A. 12) to (A.29) 
have units of slug-ft and elements (A.30) to (A.38) have units of slug-ft 2 . 
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Appendix B - Dynamic Response Comparison To A Rudder Doublet 


Nonlinear simulation results are compared with results from Boeing Commercial Aircraft Company 
and NASA Langley Research Center for the flight condition defined below. A 10 degree doublet is 
applied at all three rudder segments. The simulation was run at 40 Hz using an Euler integration 
algorithm. 


« = 1036.57435 

p = o. 

h = 20000. 0 = -1.32394e-23 

v =0. 

<7=0. 

A =0. 

6 = 0.0261307 

w = 27.092043 

r =0. 

T = 0. 

iff = 5.58502e-06 

DLEl = 10. 

DTE1 = 3. 

DTE5 = 0. 

PLA1 CP = 68.4916 

DLE2 = 0. 

DTE2 = 3. 

DTE6 = 0. 

PLA2 CP = 68.4916 

DLE3 = 0. 

DTE3 = 0. 

DTE7 = 3. 

PLA3 CP = 68.4916 

DLE4 = 10. 

DTE4 = 0. 

DTE8 = 3. 

PLA4.CP = 68.4916 

DRUD1 =0. 

DVF = 0. 

DSTAB = 

-0.426949 

DRUD2 = 0. 

DGEAR = 0. 

DELEV = 

-0.853989 


DRUD3 = 0. 

GW = 279080. 
XCGJPMAC = 54.602 
YCG_BL_MD = 0. 
ZCG_WL_MD = 201.5 

INIT_FLAG = 0 
RIGID_FLAG = 0 
ROTATION_FLAG = 0 
DRYDEN_FLAG = 0 
TLHvfiT_FLAG = 1 
AUTOCG_FLAG = 1 
AUTOWT_FLAG = 1 
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